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On the Practical use of Variable Elimination in Constraint 
Optimization Problems: 'Still-life' as a Case Study 



Variable elimination is a general technique for constraint processing. It is often dis- 
carded because of its high space complexity. However, it can be extremely useful when 
combined with other techniques. In this paper we study the applicability of variable elim- 
ination to the challenging problem of finding still-lifes. We illustrate several alternatives: 
variable elimination as a stand-alone algorithm, interleaved with search, and as a source of 
good quality lower bounds. We show that these techniques are the best known option both 
theoretically and empirically. In our experiments we have been able to solve the n = 20 
instance, which is far beyond reach with alternative approaches. 

1. Introduction 

Many problems arising in domains such as resource allocation (Gabon, de Givry, Lobjois, 
Schiex, & Warners, 1999), combinatorial auctions (Sandholm, 1999), bioinformatics and 
probabilistic reasoning (Pearl, 1988) can be naturally modeled as constraint satisfaction 
and optimization problems. The two main solving schemas are search and inference. Search 
algorithms constitute the usual solving approach. They transform a problem into a set of 
subproblems by selecting one variable and instantiating it with its different alternatives. 
Subproblems are solved applying recursively the same transformation rule. The recursion 
defines a search tree that is normally traversed in a depth-first manner, which has the 
benefit of requiring only polynomial space. The practical efficiency of search algorithms 
greatly depends on their ability to detect and prune redundant subtrees. In the worst-case, 
search algorithms need to explore the whole search tree. Nevertheless, pruning techniques 
make them much more effective. 

Inference algorithms (also known as decomposition methods) solve a problem by a se- 
quence of transformations that reduce the problem size, while preserving its optimal cost. A 
well known example is bucket elimination (BE, also known as variable elimination) (Bertele 
&: Brioschi, 1972; Dechter, 1999). The algorithm proceeds by selecting one variable at a 
time and replacing it by a new constraint which summarizes the effect of the chosen vari- 
able. The main drawback of BE is that new constraints may have large arities which require 
exponentially time and space to process and store. However, a nice property of BE is that 
its worst-case time and space complexities can be tightly bounded by a structural parame- 
ter called induced width. The exponential space complexity limits severely the algorithm's 
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practical usefulness. Thus, in the constraint satisfaction community variable elimination is 
often disregarded. 

In this paper we consider the challenging problem of finding still-lifes which are stable 
patterns of maximum density in the game of life. This academic problem has been recently 
included in the CSPlib repository^ and a dedicated web page^ has been set to maintain 
up-to-date results. In Bosch and Trick (2002), the still- life problem is solved using two 
different approaches: integer programming and constraint programming, both of them based 
on search. None of them could solve up to the n = 8 problem within reasonable time. Their 
best results were obtained with a hybrid approach which combines the two techniques and 
exploits the problem symmetries in order to reduce the search space. With their algorithm, 
they solved the n = 15 case in about 8 days of cpu. Smith (2002) proposed an interesting 
alternative using pure constraint programming techniques, and solving the problem in its 
dual form. In her work. Smith could not improve the n = 15 limit. Although not explicitly 
mentioned, these two works use algorithms with worst-case time complexity 0(2^" )). 

In this paper we show the usefulness of variable elimination techniques. First we apply 
plain BE. Against what could be expected, we observe that BE is competitive with state- 
of-the-art alternatives. Next, we introduce a more sophisticated algorithm that combines 
search and variable elimination (following the ideas of Larrosa & Dechter, 2003) and uses 
a lower bound based on mini-buckets (following the ideas of Kask &; Dechter, 2001). With 
our algorithm, we solve in one minute the n = 15 instance. We have been able to solve up to 
the n = 20 instance, which was far beyond reach with previous techniques. For readability 
reasons, we only describe the main ideas and omit algorithmic details.'^ 

The structure of the paper is the following: In the next Section we give some preliminary 
definitions. In Section 3 we solve the problem with plain BE. In Section 4 we introduce the 
hybrid algorithm with which we obtained the results reported in Section 5. In Section 6 we 
discuss how the ideas explored in this article can be extended to other domains. Besides, 
we report additional experimental results. Finally, Section 7 gives some conclusions and 
lines of future work. 

2. Preliminaries 

In this Section we first define the still-life problem. Next, we define the weighted CSP 
framework and formulate the still-life as a weighted CSP. Finally, we review the main 
solving techniques for weighted CSPS. 

2.1 Life and Still-Life 

The game of life (Gardner, 1970) is played over an infinite checkerboard, where each square 
is called a cell. Each cell has eight neighbors: the eight cells that share one or two corners 
with it. The only player places checkers on some cells. If there is a checker on it, the cell 
is alive, else it is dead. The state of the board evolves iteratively according to the following 
three rules: (1) if a cell has exactly two living neighbors then its state remains the same 

1. www.csplib.org 

2. www. ai . sri . com/" nysmith/lif e 

3. The interested reader can find an extended version, along with the source code of our implementation 
in www. Isi .upc . edu/" larrosa/publications 
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Figure 1: A: A 3 x 3 still-life. B: constraint graph of a simple WCSP instance with four 
variables and three cost functions. C: the constraint graph after assigning variable 
X4. D: the constraint graph after clustering variables and X4. E: the constraint 
graph after eliminating variable X4. 



in the next iteration. (2) if a cell has exactly three living neighbors then it is alive in the 
next iteration and (3) if a cell has fewer than two or more than three living neighbors, 
then it is dead in the next iteration. Although defined in terms of extremely simple rules, 
the game of life has proven mathematically rich and it has attracted the interest of both 
mathematicians and computer scientists. 

The still-life problem SL{n) consist on finding anxn stable pattern of maximum density 
in the game of life. All cells outside the pattern are assumed to be dead. Considering the 
rules of the game, it is easy to see that each cell must satisfy the following three 

conditions: (1) if the cell is alive, it must have exactly two or three living neighbors, (2) if 
the cell is dead, it must not have three living neighbors, and (3) if the cell is at the grid 
boundary (i.e, « = 1 or i = n or j = 1 or j = n), it cannot be part of a sequence of three 
consecutive living cells along the boundary. The last condition is needed because three 
consecutive living cells at a boundary would produce living cells outside the grid. 

Example 1 Figure l.A shows a solution to SL{3). It is easy to verify that all its cells 
satisfy the previous conditions, hence it is stable. The pattern is optimal because it has 6 
living cells and no 3 x 3 stable pattern with more that 6 living cells exists. 

2.2 Weighted CSP 

A weighted constraint satisfaction problem (WCSP) (Bistarelli, Montanari, &: Rossi, 1997) 
is defined by a tuple {X,D,T), where X = {xi, . . . ,Xn} is a set of variables taking values 
from their finite domains Di ^ D. is a set of weighted constraints [i.e., cost functions). 
Each / G is defined over a subset of variables, var{f), called its scope. The objective 
function is the sum of all functions in !F, 

and the goal is to find the instantiation of variables that minimizes the objective function. 

Example 2 Consider a WCSP with four variables X = {xi}f^i with domains Di = {0, 1} 
and three cost functions: fi{xi^Xji) = xi +2:4, /2(a;2,a;3) = X2X^ and f3{x2,X4) = X2 + X4. 



423 



Larrosa, Morancho & Niso 



The objective function is F{xi,X2,X3,X4) = xi + X4 + X2X3 + X2 + X4. Clearly, the optimal 
cost is 0, which is obtained with every variable taking value 0. 

Constraints can be given explicitly by means of tables, or implicitly as mathematical 
expressions or computing procedures. Infeasible partial assignments are specified by con- 
straints that assign cost 00 to them. The assignment of value a to variable Xi is noted 
Xi = a. A partial assignment is a tuple t = {xi^ = vi,Xi^ = V2, - • • , xi- = Vj). The extension 
of t to Xi = a is noted t ■ {xi = a). WCSPs instances are graphically depicted by means 
of their interaction or constraint graph, which has one node per variable and one edge con- 
necting any two nodes that appear in the same scope of some cost function. For instance, 
Figure l.B shows the constraint graph of the problem in the previous example. 

2.3 Overview of Some Solving Techniques 

In this Subsection we review some solving techniques widely used when reasoning with 
constraints. 

2.3.1 SEARCH 

WCSPs are typically solved with depth-first search. Search algorithms can be defined in 
terms of instantiating functions. 

Definition 1 Let P = {X,D,J^) a WCSP instance, f a function in T, Xi a variable in 
var{f), and v a value in Di. Instantiating / with Xi = v is a new function with scope 
var[f) — {xi} which returns for each tuple t, f{t ■ [xi = v)). Instantiating P with Xi= v is 
a new problem P \x.-^= [X — {xi}, D — {Di}, T'), where T' is obtained by instantiating all 
the functions in T that mention Xi with Xi = v. 

For instance, instantiating the problem of Example 2 with X4 = 1, produces a new 
problem with three variables {xi}f^i and three cost functions: fi{xi,X4 = 1) = a;i + 1, 
f2{x2,X3) = X2X3 and f3{x2,X4 = 1) = X2 + 1. Figure l.C shows the corresponding 
constraint graph, obtained from the original graph by removing the instantiated variable X/i 
and all adjacent edges. Observe that the new graph depends on the instantiated variable, 
but does not depend on the value assigned to it. 

Search algorithms transform the current problem P into a set of subproblems. Usually 
it is done by selecting one variable Xi which is instantiated with its different domain values 
(P |xi=t;i5-^ U,=t;2? ■ ■ ■ ! \x,=vd)- This transformation is called branching. In each subprob- 
lem the same process is recursively applied, which defines a tree of subproblems. Search 
algorithms expand subproblems until a trivial case is achieved: there is no variable left, or 
a pruning condition is detected. In optimization problems, pruning conditions are usually 
defined in terms of lower and upper bounds. Search keeps the cost of the best solution so 
far, which is an upper bound of the optimal cost. At each node, a lower bound of the best 
cost obtainable underneath is computed. If the lower bound is greater than or equal to the 
upper bound, it is safe to backtrack. 

The size of the search tree is 0{d" ) (being d the size of the largest domain) which bounds 
the time complexity. If the tree is traversed depth-first, the space complexity is polynomial. 
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2.3.2 CLUSTERING 

A well-known technique for constraint processing is clustering (Dechter & Pearl, 1989). It 

merges several variables into one meta-variable, while preserving the problem semantics. 
Clustering variables Xi and Xj produces meta- variable x^, whose domain is Di x Dj. Cost 
functions must be accordingly clustered. For instance, in the problem of Example 2, cluster- 
ing variables 2:3 and X4 produces variable Xc with domain Dc = {(0, 0), (0, 1), (1, 0), (1, 1)}. 
Cost fimctions /2 and /a are clustered into fc{x2-,Xc) = /2 + /a- With the new variable 
notation /c = 2;2a;c[l] + X2 + Xc\^], where Xc[i\ denotes the 'i-th component of Xc- Function 
/i needs to be reformulated as fi{xi,Xc) = xi +xj(2]. The constraint graph of the resulting 
problem is obtained by merging the clustered variables and connecting the meta-node with 
all nodes that were adjacent to some of the clustered variables. Figure l.D shows the con- 
straint graph after the clustering of x^ and X4. The typical use of clustering is to transform 
a cyclic constraint graph into an acyclic one, which can be solved efficiently thereafter. 

2.3.3 VARIABLE elimination 

Variable elimination is based on the following two operations. 

Definition 2 The sum of two functions f and g, noted [f + g), is a new function with 
scope var{f) U var{g) which returns for each tuple the sum of costs of f and g, 

{f + 9){t)=f{t)+9{t) 

Definition 3 The elimination of variable xi from f, noted f ij. xi, is a new function with 
scope var{f) — {xi} which returns for each tuple t the cost of the best extension oft to Xi, 

{f^Xi){t) = min{/(i-(rEi = a))} 

Observe that when / is a unary function (i.e., arity one), eliminating the only variable 
in its scope produces a constant. 

Definition 4 Let P = {X,D,T) be a WCSP instance. Let Xi E X be an arbitrary variable 
and let Bi be the set of all cost functions having Xi in their scope (Bi is called the bucket of 
Xi). We define gi as 

9i = (Y. f") 

The elimination of Xi transforms P into a new problem P JJ-xi= {X — {xi\,D — {-Dj}, — 
Bi) U {ffj}}. In words, P is obtained by replacing Xi and all the functions in its bucket 
by gi. 

P and P ij-xi have the same optimal cost because, by construction, gi compensates the 
absence of Xi. The constraint graph of P ij-xi is obtained by forming a clique with all the 
nodes adjacent to node Xi and then removing Xi and all its adjacent edges. For example, 
eliminating X4 in the problem of Example 2 produces a new problem with three variables 
{xi}^^^ and two cost functions: /2 and g^. The scope of 34 is {xi,X2} and it is defined as. 
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94 = (/i + /s) JJ- = (2:1 + X4 + X2 + X4) i}. X4 = xi + X2. Figure l.D shows the constraint 
graph after the elimination. 

In the previous example, the new function could be expressed as a mathematical ex- 
pression. Unfortunately, in general, the result of summing functions or eliminating variables 
cannot be expressed intensionally, and new cost functions must be stored extensionally in 
tables. Consequently, the space complexity of computing P ij-^^ is proportional to the num- 
ber of entries of gi, which is: &{YlxjGvar{gi) l-^jD- Since Xj G var{gi) iS Xj is adjacent to Xi 
in the constraint graph, the previous expression can be rewritten as @{YlxjeN{i,Gp) 
where Gp is the constraint graph of P and N{i, Gp) is the set of neighbors of Xi in Gp. 
The time complexity of computing P JJ-^;. is its space complexity multiplied by the cost of 
computing each entry of gi. 

Bucket elimination (BE) works in two phases. In the first phase, it eliminates variables 
one at a time in reverse order. In the elimination of Xi, the new gi function is computed 
and added to the corresponding bucket. The elimination of xi produces an empty-scope 
function (i.e., a constant) which is the optimal cost of the problem. In the second phase, BE 
considers variables in increasing order and generates the optimal assignment of variables. 
The time and space complexity of BE is exponential on a structural parameter from the 
constraint graph, called induced width, which captures the maximum arity among all the 
gi functions. Without any additional overhead BE can also compute the number of optimal 
solutions (see Dechter, 1999, for details). 

2.3.4 SUPER-BUCKETS 

In some cases, it may be convenient to eliminate a set of variables simultaneously (Dechter 
& Fatah, 2001). The elimination of the set of variables Y is performed by collecting in By 
the set of functions mentioning at least one variable of Y. Variables in Y and functions in 
By are replaced by a new function gy defined as. 

The set By is called a super-bucket. Note that the elimination of Y can be seen as the 
clustering of its variables into a meta-variable xy followed by its elimination. 

2.3.5 MINI-BUCKETS 

When the space complexity of BE is too high, an approximation, called mini buckets 
(Dechter Sz Rish, 2003), can be used. Consider the elimination of Sj, with its associated 
bucket Bi = {fi^, . . . , /jj.}. BE would compute, 

9i = {J2 

The time and space complexity of this computation depends on the arity of gi. If it is beyond 
our available resources, we can partition bucket Bi into so-called mini-buckets Bi^^ . . . , Bi^, 
where the number of variables in the scopes of each mini-bucket is bounded by a parameter. 
Then we can compute. 
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Figure 2: A constraint graph and its evolution over a sequence of variable eliminations and 
instantiations. 

where each gi^ has a bounded arity. Since, 

feBi j=i feBi. 

the elimination of variables using mini-buckets yields a lower bound of the actual optimal 
cost. 

2.3.6 COMBINING SEARCH AND VARIABLE ELIMINATION 

When plain BE is too costly in space, we can combine it with search (Larrosa & Dechter, 
2003). Consider a WCSP whose constraint graph is depicted in Figure 2. A. Suppose that 
we want to eliminate a variable but we do not want to compute and store constraints with 
arity higher than two. Then we can only take into consideration variables connected to at 
most two variables. In the example, variable is the only one that can be selected. Its 
elimination transforms the problem into another one whose constraint graph is depicted in 
Figure 2.B. Now xq has its degree decreased to two, so it can also be eliminated. The 
new constraint graph is depicted in Figure 2.C. At this point, every variable has degree 
greater than two, so we switch to a search schema which selects a variable, say 2:3 , branches 
over its values and produces a set of subproblems, one for each value in its domain. All of 
them have the same constraint graph, depicted in Figure 2.D. For each subproblem, it is 
possible to eliminate variable xg and X4. After their elimination it is possible to eliminate 
X2 and xg, and subsequently X5 and xi. Eliminations after branching have to be done at 
every subproblem since the new constraints with which the eliminated variables are replaced 
differ from one subproblem to another. In the example, only one branching has been made. 
Therefore, the elimination of variables has reduced the search tree size from to d, where 
d is the size of the domains. In the example, we bounded the arity of the new constraints 
to two, but it can be generalized to an arbitrary value. 

3. Solving Still-life with Variable Elimination 

SL{n) can be easily formulated as a WCSP. The most natural formulation associates one 
variable xij with each cell («, j). Each variable has two domain values. If xij = the cell is 
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Figure 3: A: Structure of the constraint graph oi SL{n). The node in the center, associated 
to cell is linked to all cells it interacts with. The shadowed area indicates 

the scope of /y . B (left): Constraint graph of 5L(6) after clustering cells into 
row variables. B (from left to right: Evolution of the constraint graph during the 
execution of BE. 



dead, if Xij = 1 it is alive. There is a cost function for each variable xij. The scope of 
fij is Xij and all its neighbors. It evaluates the stability of Xij: if Xij is unstable given its 
neighbors, fij returns oo; else fij returns 1 — Xij."^ The objective function to be minimized 
is, 

n n 

If the instantiation X represents an unstable pattern, F{X) returns oo; else it returns the 
number of dead cells, fij can be stored as a table with 2^ entries and evaluated in constant 
time. 

Figure 3.^4 illustrates the structure of the constraint graph of SL{n). The picture shows 
an arbitrary node Xij linked to all the nodes it interacts with. For instance, there is an edge 
between Xij and Xij^i because Xij^i is a neighbor of Xij in the grid and, consequently, 
both variables are in the scope of fij. There is an edge between Xij and Xi-ij-2 because 
both cells are neighbors of Xi^ij^i in the grid and, therefore, both appear in the scope of 
fi-ij-i. The shadowed area represents the scope of fij (namely, Xij and all its neighbors). 
The complete graph is obtained by extending this connectivity pattern to all nodes in the 
graph. 

For the sake of clarity, we use an equivalent but more compact SL{n) formulation 
that makes BE easier to describe and implement: we cluster all variables of each row 
into a single meta-variable. Thus, Xi denotes the state of cells in the i-th row (namely, 
Xi = {xii,Xi2, ■ ■ ■ ,Xin) with Xij G {0, 1}). Accordingly, it takes values over the sequences of 
n bits or, equivalently, over the natural numbers in the interval [0..2" — 1]. Cost functions 
are accordingly clustered: there is a cost function fi associated with each row i, defined as, 

n 

fi = E 

i=i 

4. Recall that, as a WCSP, the task is to minimize the number of dead cells. Therefore, we give cost 1 to 
dead cells and cost to living cells. 
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For internal rows, tliG scope of fi is ^Xi—i^Xi^Xi^i }. The cost function of the top row, /i, 
has scope {xi,X2}. The cost function of the bottom row, /„, has scope {.t„_i,.t„}. If there 
is some unstable cell in Xi, fi{xi-i,Xi,Xi+i) = oo. Else, it returns the number of dead cells 
in Xi. Evaluating fi is 0(n) because all the bits of the arguments need to be checked. The 
new, equivalent, objective function is, 

n 

Figure 3.S (left) shows the constraint graph of SL{Q) with this formulation. An arbitrary 
variable Xi is connected with the two variables above and the two variables below. The 
sequential structure of the constraint graph makes BE very intuitive. It eliminates variables 
in decreasing orders. The elimination of Xi produces a new function gi = {fi-i + gi+i) i Xi 
with scope {xi^2,Xi^i}. Figure 3.B (from left to right) shows the evolution of the constraint 
graph along the elimination of its variables. Formally, BE applies a recursion that transforms 
subproblem P into P JJ-^; . , where Xi is the variable in P with the highest index. It satisfies 
the following property. 

Property 1 Let gi be the function added by BE to replace Xi. Then gi{a,b) is the cost of 
the best extension of {xi-2 = a, Xi-i = b) to the eliminated variables {xi, . . . ,Xn)- Formally, 

gi{a,b) = min {fi^i{a,b,Vi) + fi{b,Vi,Vi+i) + 

VieDi,...,VneDn 

+ fi+l{vi,Vi+i,Vi+2) + ■ ■ ■ 

+fn-l{Vn-2,Vn-l,Vn) + fn{Vn-l,Vn)} 

If gi{a,b) = oo, it means that the pattern a,b cannot be extended to the inferior rows 
with a stable pattern. If gi{a, b) = k (with k / oo), it means that a, b can be extended and 
the optimal extension has k dead cells from Xi-i to Xn- 

The space complexity of BE ©(n x 2^"), due to the space required to store n functions 
gi extensionally (2" x 2" entries each). Regarding time, computing each entry of gi has 
cost 6(n X 2") (finding the minimum of 2" alternatives, the computation of each one is 
6(n)). Since each gi has 2^" entries, the total time complexity is 0(n^ x 2^"). Observe that 
solving SL{n) with BE is an exponential improvement over search algorithms, which have 
time complexity 0(2" ). 

Table 4 reports some empirical results. They were obtained with a 2 Ghz Pentium IV 
machine with 2 Gb of memory. The first columns reports the problem size, the second 
reports the optimal cost as the number of dead cells (in parenthesis, the number of living 
cells), the third column reports the number of optimal solutions. We count as different 
two solutions even if one can be transformed to the other through a problem symmetry. 
The fourth column reports the CPU time of BE in seconds. The fifth, sixth and seventh 
columns report the results obtained with the three approaches tried by Bosch and Trick 
(2002):^ constraint programming (CP), integer programming (IP), and a more sophisticated 
algorithm (CP/IP) which combines CP and IP, and exploits the problem symmetries. 

5. The corresponding OPL code is available at http://mat.gsia.cmu.edu/LIFE. 
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Figure 4: Experimental results of four different algorithms on the still-life problem. Times 
are in seconds. 



It can be observed that BE clearly outperforms CP and IP by orders of magnitude. 
The n = 14 case is the largest instance that we could solve due to exhausting the available 
space. Comparing BE with CP/IP, we observe that there is no clear winner. An additional 
observation is that BE scales up very regularly, each execution requiring roughly eight times 
more time and four times more space than the previous, which is in clear accordance with 
the algorithm complexity. 

4. Combining Search and Variable Elimination 

One way to overcome the high space complexity of BE is to combine search and variable 
elimination in a hybrid approach HYB (Larrosa &; Schiex, 2003). The idea is to use search 
(i.e, instantiations) in order to break the problem into independent smaller parts where 
variable elimination can be efficiently performed. 

Let us reformulate the problem in a more convenient way for the hybrid algorithm. For 
the sake of simplicity and without loss of generality consider that n is even. We cluster 
row variables into three meta- variables: xf denotes the two central cells of row i, xf and 
xf denote the f — 1 remaining cells on the right and left, respectively (see Figure 5.^). 
Consequently, x'^ takes values in the range [0..3], x^ and x^ take values in the range 
[0..2'2~^ — 1]. Cost functions are accordingly clustered, 

2 n 
fi = ^ fij ? fi = 5Z f'-j 

The new, equivalent, objective function is, 

n 
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Figure 5: Formulation of SL{n) used by the hybrid algorithm. A: Each row is clustered 
into three variables. B: Constraint graph of S'L(6). C: Constraint graph after 
the assignment of x^, x'^_i and x'^_2. D: Constraint graph after the elimination 



of and x^. 
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The scopes of internal row functions, and f-^. are a;^;^, , xf, xf^^, and 

{xf_-^,xf_i,xf ,xf,xf_^Y^^'i+i\- Top functions and have scopes {xi,xf ,X2,X2} 
and {rrf , rrf , 2:2^, a;2^}. Bottom functions and have scopes {x^_i,x^_i,x^,x^} and 
a;^_^, 2;^, 2;^}. Figure 5.S shows the corresponding constraint graph. The impor- 
tance of this formulation is that x^ and xf are independent (i.e, there is no edge in the 
constraint graph connecting left and right variables). 

The hybrid algorithm HYB searches over the central variables and eliminates the lateral 
variables. Variables are considered in decreasing order of their index. Thus, the algorithm 
starts instantiating x^, x^_i and x^_2, which produces a subproblem with the constraint 
graph shown in Figure 5.C Observe that variable x^ (respectively, 2;^) is only connected 
with variables x^_i and x^_2 (respectively, x^_i and a;^_2). Then it is eliminated producing 
a new function with scope {x^-2T^n-i} (respectively, with scope {x!^-2j^n-i})- 
Figure 5.D shows the resulting constraint graph. Lateral variables have domains of size 
2^2 ^ . Hence, their elimination is space 6(2") and time e (2^1 ). It is important to note that 
these eliminations are subject to the current assignment of a;^, x^_i and x^_2. Therefore, 
they have to be recomputed when their value change. After the elimination of x^ and 
x^, the algorithm would assign variable x^_^ which will make possible the elimination of 
x^_i and x^_i, and so on. At an arbitrary level of search, the algorithm assigns xf, which 
makes xf_^2 ^^i^l 2;,^2 independent of the central columns and only related to their two 
variables above. Then, it eliminates them by replacing the variables by functions 5/^2 ^.nd 
9i+2 with scopes {xf,xf_^^} and {a;|*, a;^^^}, respectively. Formally, HYB applies a recursion 
that transforms subproblem P into 4 simpler subproblems {{{P Lc^„) ij-^L ) ij- r }l^o- It 
satisfies the following property. 



Property 2 Let gf he the function computed by HYB used to replace variable xf. Then 
gf'{a,b) is the cost of the best extension of {xf_2 = a, xf_^ = b) to eliminated variables 
(a;f , . . . , a;^), conditioned to the current assignment. Similarly, for the right side, gf{a,b) 
is the cost of the best extension of {xf_2 = a, a;fi^ = b) to eliminated variables {xf, . . . , a;^), 
conditioned to the current assignment. 



A consequence of the previous Property is that the minimum 52^2 («i b) among all com- 
binations of a and 6 is a lower bound of the best cost that can be obtained in the left 
part of the grid if we continue the current line of search. Therefore, m.ma^b{9i+2i^^^)} + 
m.inafi{gi^2i^^^)} ^ valid lower bound of the current node and can be used for pruning 
purposes. 

The space complexity of the algorithm is 0(n x 2"), due to the g[ and gf^ functions 
which need to be explicitly stored. The time complexity is 0(n x 2^-^"), because 0(4") 
nodes may be visited (n variables with domains of size 4) and the cost of processing each 
node is 6(n X 2^2') due to the variable eliminations. 

Thus, comparing with BE, the time complexity increases from ©(n^ x2^") to 0(nx2^'^"). 
This is the prize HYB pays for the space decrement from 0(n x 2^") to 0(n x 2"). 
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4.1 Refining the Lower Bound 

It is well-known that the average-case efficiency of search algorithms depends greatly on 
the lower bound that they use. Our algorithm is using a poor lower bound based on the g-" 
and g-^ functions, only. 

Kask and Dechter (2001) proposed a general method to incorporate information from 
yet-unprocessed variables into the lower bound. Roughly, the idea is to run mini buckets 
(MB) prior search and save intermediate functions for future use. MB is executed using the 
reverse order in which search will instantiate the variables. When the execution of MB is 
completed, the search algorithm is executed. At each node, it uses mini-bucket functions 
as compiled look-ahead information. In this Subsection, we show how we have adapted this 
idea to SL{n) and how we have integrated it into HYB. 

Consider SL{n) formulated in terms of left, central and right variables [xf^xf^xf-). 
The exact elimination of the first row variables [x^^x^ ^x^) can be done using super-bucket 
Bi = {fit fi', f2 1 f2'} ^^'^ computing the function, 

hi = {ft + ft + fi + fi) ^ 4} 

The scope oi hi is {x2,X2 ,X2,x]^,x^ ,xf}. Using the mini-buckets idea, we partition the 
bucket into = {ft if 2} and = {ft^,fi}- Then, we approximate hi by two smaller 
functions hf and hf, 

hf = {ft + fi)i^{^t:^?} 
hf = {fl^ + f2^)i^{^?:xf} 

The scopes of /if and /if are {x2,X2 , x^,x^} and {x2 ,X2, x^ , x^}, respectively. The same 
idea is repeated row by row in increasing order. In general, processing row i, yields two 
functions, 

hf = {hf_,+ft+i)i^{xf,xf} 
hf = ihf_, + ftl,)i^{xf,xt} 

The scopes of hf and hf are {xf_^^, xf_^^, xt_^2^ xf^2} and {xf_^^,xt_^-^,xf_^2^xt^2}^ respec- 
tively. By construction, ht{a,a' ,b,b') contains the cost of the best extension of a, a', 6, 6' 
to processed variables xt,xf',..., x^^xt considering left functions only. We have the same 
property for ht{a' ,a,b' ,b) and right functions. 

The complexity of MB is space B(nx2") and time 0(n^ x2^'^"). Since these complexities 
are smaller than the complexity of HYB, running this pre-process does not affect its overall 
complexity. 

After MB is executed, HYB can use the information recorded in the /i,f and hf functions. 
Consider an arbitrary node in which HYB assigns xf and eliminates 2:^2 and xl^2- Let a 
and b be domain values of variables and Prom Property 2 we have that 5/^2 (^^5 ^) 

contains the best extension of a, b that can be attained in the left part of rows i + 1 
to n as long as the current assignment X*^ is maintained. Additionally, we have that 
ht_i{a,xf ,b,xf_^i) contains the best extension of a, 6 that can be attained in the left part 
of rows « to 1. Therefore, gt^2 

(a, 6) + /if_i(a J is a lower bound for a,b and X 

of the left part of the grid. Consequently, 

■"^°a,he[o..2t-i-i]{5^H2(«'^) +ht-i{a,xf,b,xf^^)} 
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is a lower bound of the left part of the grid for the current assignment. With the same 
reasoning on the right part we have that, 

■"^^a,6e [0..2 1 - 1 - 1] {9^+2 {a,b) + hf_^{a,xf ,b,xf^^)} + 
+'"^°a,fee[o..2t-i-i]{f 1^2(0: b) + h^_i{xf, a, xf^^,b)} 
is a lower bound of the current assignment. 

4.2 Refining the Upper Bound 

The efficiency of the algorithm also depends on the initial value of the upper bound. A 
good upper bound facilitates pruning earlier in the search tree. Bosch and Trick (2002) 
suggested to modify SL{n) by adding the additional constraint of considering symmetric 
patterns, only. Since the space of solutions becomes considerably smaller, the problem is 
presumably simpler. Clearly, the cost of an optimal symmetric stable pattern is an upper 
bound of the optimal cost of SL{n). It has been observed that such upper bounds are very 
tight. 

Since the motivation of our work is to use variable elimination techniques, we have 
considered still-lifes which are symmetric over a vertical reflection, because they can be 
efficiently solved using BE. The symmetric still-life problem SSL{n) consists on finding a 
n X n stable pattern of maximum density in the game of life subject to a vertical reflection 
symmetry (namely, the state of cells and {i,n — j + 1) must be the same.^ 

Adapting BE to solve SSL{n) is extremely simple: we only need to remove symmetrical 
values from the domains. Let us assume that n is an even number (the odd case is similar). 
We represent a symmetric sequences of bits of length n by considering the left side of the 
sequence (i.e, the first n/2 bits). The right part is implicit in the left part. Thus, we 
represent symmetrical sequences of n bits as integers in the interval [0..2'2 — 1]. Reversing a 
sequence of bits a is noted a. Hence, if a is a sequence of n/2 bits, a ■ a is the corresponding 
symmetrical sequence of n bits. 

The complexity of BE, when applied to SSL{n) is time 6(n^ x2^-^") and space 6(nx2"). 
Therefore, executing it prior HYB and setting the upper bound with its optimal cost does 
not affect the overall complexity of the hybrid. 

4.3 Further Exploitation of Symmetries 

SL{n) is a highly symmetric problem. For any stable pattern, it is possible to create an 
equivalent pattern by: (i) rotating the board by 90, 180 or 270 degrees, {ii) reflecting 
the board horizontally, vertically or along one diagonal or {iii) doing any combination of 
rotations and reflections. 

Symmetries can be exploited at very different algorithmic levels. In general, we can 
save any computation whose outcome is equivalent to a previous computation due to a 
symmetry if we have kept its outcome. For instance, in MB it is not necessary to compute 
hf{a\a,b' ,b) because it is equal to hf{a,a',b,b') due to the vertical reflection symmetry. 
Another example occurs in HYB. Let x^ = Vn,x^_i = Vn-i, ■ ■ ■ ,xf = Vi he the current 

6. Unlike Smith's (2002) work we cannot easily exploit a larger variety of symmetries such as rotations and 
diagonal reflections. 



434 



On the practical use of variable elimination 



n 


opt 


opt-SSL 


CP/IP 


BE 


HYB 


HYB no LB 


HYB no UB 


13 


79(90) 


79 


12050 


13788 


2 


2750 


2 


14 


92(104) 


92 


5 X 10° 


10^ 


2 


7400 


3 


15 


106(119) 


106 


7 X 10^ 


* 


58 


2 X 10^ 


61 


16 


120(136) 


120 


* 


* 


7 


6 X 10^ 


49 


17 


137(152) 


137 


* 


* 


1091 


* 


2612 


18 


153(171) 


154 


* 


* 


2029 


* 


2311 


19 


171(190) 


172 


* 


* 


56027 


* 


56865 


20 


190(210) 


192 


* 


* 


2 X 10^ 


* 


2 X 10^ 


22 


? 


232 


* 


* 


* 


* 


* 


24 


? 


276 


* 


* 


* 


* 


* 


26 


? 


326 


* 


* 


* 


* 


* 


28 


? 


378 


* 


* 


* 


* 


* 



Figure 6: Experimental results of three different algorithms on the still-life problem. Times 
are in seconds. 



assignment. The reversed assignment = v^, x^-i = Wjj-i, . . . , ai^" = i)J is equivalent due 
to the vertical reflection symmetry. Thus, if it has already been considered, the algorithm 
can backtrack. Our implementation uses these tricks and some others which we do not 
report because it would require a much lower level description of the algorithms. 

5. Experimental Results 

Figure 6 shows the empirical performance of our hybrid algorithm. The first column contains 
the problem size. The second column contains the optimal value as the number of dead 
cells (in parenthesis the corresponding number of living cells). The third column contains 
the optimal value of the symmetrical problem SSL{n), obtained by executing BE. It can 
be observed that SSL{n) provides very tight upper bounds to SL{n). The fourth column 
reports the time obtained with the CP/IP algorithm (Bosch & Trick, 2002). The fifth 
column reports times obtained with BE. The sixth column contains times obtained with 
our hybrid algorithm HYB. As it can be seen, the performance of HYB is spectacular. The 
n = 14 and n = 15 instances, which require several days of CPU, are solved by HYB in a few 
seconds. Instances up to n = 18 are solved in less than one hour. The largest instance that 
we can solve is n = 20, which requires about two days of CPU (Figure 7 shows the optimal 
n = 19 and n = 20 still- lifes) . Regarding space, our computer can handle executions of 
HYB up to n = 22. However, neither the n = 21 nor the n = 22 instance could be solved 
within a week of CPU. It may seem that solving the n = 20 instance is a petty progress with 
respect previous results on the problem. This is clearly not the case. The search space of 
the n = 15 and n = 20 instances have size 2^^ = 2^^^ and 2^° = 2^°°, respectively. Thus, 
we have been able to solve a problem with a search space 2^^^ times larger than before. 
Since BE scales up very regularly, we can accurately predict that it would require 4000 Gb 
of memory and about 7 centuries to solve the n = 20 instance. 
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Figure 7: Maximum density still-lifes for n = 19 and n = 20. 



Since HYB combines several techniques, it is interesting to assess the impact of each 
one. The seventh column reports times obtained with HYB without using mini-buckets 
information in the lower bound. As can be seen, the algorithm is still better than plain BE, 
but it performance is dramatically affected. The information gathered during the preprocess 
improves the quality of the lower bound and anticipates pruning. Finally, the eighth column 
reports times obtained with HYB without having the upper bound initialized to SSL{n). 
In this case we see that the importance of this technique is quite limited. The reason is 
that HYB, even with a bad initial upper bound, finds the optimum very rapidly and, after 
that moment, the quality of the initial upper bound becomes irrelevant. 

6. Extension to Other Domains 

The SL(n) problem has a very well defined structure, and the hybrid algorithm that we 
have proposed makes an ad hoc exploitation of it. It is easy to find the right variables to 
instantiate and eliminate. It is also easy to find a variable order for which mini buckets 
produces good quality lower bounds. A natural question is whether it is possible to apply 
similar ideas to not so well structured problems. The answer is that it is often possible, 
although we need to rely on more naive and consequently less efficient exploitation of the 
problems' structure. In this Section we support our claim by reporting additional exper- 
imental results on different benchmarks. In particular, we consider spotS and DIMACS 
instances. Spot5 instances are optimization problems taken from the scheduling of an earth 
observation satellite (Bensana, Lemaitre, Sz Verfaillie, 1999). The DIMACS benchmark con- 
tains SAT instances from several domain. Since we are concerned with optimization tasks, 
we have selected some unsatisfiable instances and solved the Max-SAT task (i.e, given an 
unsatisfiable SAT instance, find the maximum number of clauses that can be simultaneously 
satisfied), which can be modeled as a WCSP (de Givry, Larrosa, Meseguer, & Schiex, 2003). 
We consider aim instances (artificially generated random 3-SAT), pret (graph coloring), ssa 
and 6/ (circuit fault analysis). 

Figure 8 shows the constraint graph of one instance of each domain, as visualized by 
LEDA graph editor. It can be observed that these graphs do not have an obvious pattern 
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Figure 8: Constraint graph of four WCSP instances. Prom the top-left corner, clockwise, 
aim-lOO-l-6-no-l, pret60-25, ssa0432-003 and Spot5-404. 



to be exploited. Thus, we have to use variable elimination techniques in a more naive way. 
We solve the problems with the generic WCSP solver TOOLBAR^ (TB). It performs a depth- 
first branch-and-bound search and it is enhanced with general-purpose dynamic variable and 
value ordering heuristics. We modified TOOLBAR to combine search and variable elimination 
as follows: at an arbitrary subproblem, every variable with degree less than 3 is eliminated. 
Only when all the variables have degree larger than or equal to 3, an unassigned variable is 
heuristically selected and each of its domain values are heuristically ordered and sequentially 
instantiated. The process is recursively applied to each of the subproblems. Note that 
this is a generic version of the HYB algorithm where the decision of which variables are 
instantiated and which variables are eliminated is left to a heuristic, instead of establishing 

7. Available at http://carlit.toulouse.inra.fr/cgi-bin/awki.cgi/SoftCSP. 
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it by hand. We will refer to this implementation as TB hyb- Toolbar offers a variety 
of lower bounds based on different forms of local consistency (Larrosa & Schiex, 2003). 
One of them, directional arc consistency (DAC*), is essentially equivalent to mini-buckets 
of size 2 and, therefore, similar in spirit to the lower bound computed by HYB. However, 
unlike HYB where mini-buckets are executed only once as a pre-process, toolbar executes 
DAC* at every search state, subject to the current subproblem. It has been shown by Kask 
(2000) that this approach is generally more efficient. The other main difference with respect 
HYB, is that TOOLBAR executes DAC* subject to an arbitrary variable ordering (in HYB 
a good order was identified from the problem structure). Other lower bounds available 
in TOOLBAR are node consistency (NC*) which is weaker than DAC*, and full directional 
arc consistency (FDAC*) which can be seen as a (stronger) refinement of DAC*. We have 
experimented with four algorithms: TB^"^*, TB^^^"^*, TBg^^g* and TB^^^*^*, where A^ 
denotes algorithm A with lower bound B. 

Most spots instances are too difficult for TOOLBAR. Therefore, we decreased their size 
by letting TOOLBAR make a sequence of k greedy assignments driven by its default variable 
and value ordering heuristics. The result is a subproblem with k less variables. In the 
following, Ik denotes instance / where k variables have been greedily assigned by toolbar 
with default parameters. 

Table 9 reports the result of these experiments. The first column indicates the instances 
and subsequent columns indicate the CPU time (in seconds) required by the different algo- 
rithms. A time limit of 3600 seconds was set up for each execution. It can be observed that 
TOOLBAR with the weakest lower bound (TB^^*) is usually the most inefficient alternative. 
It cannot solve any of the spotS instances and also fails with several aim and ssa instances. 
When TOOLBAR is enhanced with a mini buckets lower bound (tb^'^'"*) all spotS problems 
are solved. In the other domains, the new lower bound does not produce a significant ef- 
fect. When we further add variable elimination (TB^y^*) all the problems are solved. In 
general, there is a clear speed-up. The worst improvements are in the pret instances where 
the time is divided by a factor of 2 and the best ones are obtained in the spot5 5034o and 
ssa7552-158 instances which are solved instantly. Typical speed-ups range from 5 to 10. 
Finally, we observe that the addition of the stronger lower bound (TB^y^*^*) has a limited 
effect in these problems. Only the execution of instance ssa7552-038 is clearly accelerated. 
Therefore, from these experiments we can conclude that the main techniques that we used 
to solve the still-life problem can also be successfully applied to other domains. 

7. Conclusions 

In this paper we have studied the applicability of variable elimination to the problem of 
finding still-lifes. Finding still-lifes is a challenging problem and developing new solving 
techniques is an interesting task per se. Thus, the first contribution of this paper is the 
observation that plain variable elimination (i.e, BE) is competitive in practice and provides 
time complexity exponentially better than search-based approaches. Besides, we have de- 
veloped an algorithm with which we have been able to solve up to the n = 20 instance, 
with which we clearly improved previous results. The second contribution of the paper 
has a deeper insight. Our algorithm uses recent techniques based on variable elimination. 
Since these techniques are little known and rarely applied in the constraints community. 
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Figure 9: Experimental results in some WCSP instances with four different algorithms. 

Each column reports CPU time in seconds. Symbol - indicates that a time limit 
of 3600 seconds has been reached. 



the results presented in this paper add new evidence of their potential. We have also shown 
that variable elimination can be used beyond the academic still-life problem by providing 
experimental results in some unstructured realistic problems from different domains. 
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